Dicke Quantum Spin and Photon Glass in Optical Cavities: 
Non-equilibrium theory and experimental signatures 
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In the context of ultracold atoms in multimode optical cavities, the appearance of a quantum-critical glass 
phase of atomic spins has been predicted recently. Due to the long-range nature of the cavity-mediated interac- 
tions, but also the presence of a driving laser and dissipative processes such as cavity photon loss, the quantum 
optical realization of glassy physics has no analog in condensed matter, and could evolve into a "cavity glass 
microscope" for frustrated quantum systems out-of-equilibrium. Here we develop the non-equilibrium theory of 
the multimode Dicke model with quenched disorder and Markovian dissipation. Using a unified Keldysh path 
integral approach, we show that the defining features of a low temperature glass, representing a critical phase 
of matter with algebraically decaying temporal correlation functions, are seen to be robust against the presence 
of dissipation due to cavity loss. The universality class however is modified due to the Markovian bath. The 
presence of strong disorder leads to an enhanced equilibration of atomic and photonic degrees of freedom, in- 
cluding the emergence of a common low-frequency effective temperature. The imprint of the atomic spin glass 
physics onto a "photon glass" makes it possible to detect the glass state by standard experimental techniques of 
quantum optics. We provide an unambiguous characterization of the superradiant and glassy phases in terms of 
fluorescence spectroscopy, homodyne detection, and the temporal photon correlation function g (2) (r). 

PACS numbers: 37.30+i, 42.50.-p, 05.30.Rt, 75.10.Nr 



I. INTRODUCTION 

An emerging theme in the research on strongly correlated 
ultracold atoms is the creation of quantum soft matter phases 
ranging from nematics and smectics H |2), liquid crystals 
1151 , granular materials EH6), friction phenomena in nonlin- 
ear lattices 17][E)> to glasses [9-13 1. Realizing glasses with 
strongly interacting light-matter systems bears the promise to 
study some of the most celebrated achievements in statistical 
mechanics from a new vantage point. The Parisi solution of 
mean-field spin glasses [14|, for example, continues to trig- 
ger research more than three decades after its discovery in the 
early 1980's, and may have implications for information stor- 
age [15 1 and "frustrated" optimization algorithms |16|. The 
latter is related to the inability of a glass to find its ground 
state; a feature that makes it inherently non-equilibrium. 

Historically, quantum effects in soft matter and glasses have 
not played a prominent role because most soft materials are 
too large, too heavy and/or too hot and therefore way out- 
side the quantum regime. Spin and charge glass features have 
however been invoked in some electronic quantum materials 
lfl4l[T7ll . mainly due to RKKY-type interactions or randomly 
distributed impurities providing a random potential for the 
electrons. However, here the glassy mechanisms occur of- 
ten in combination with other more dominant (Coulomb) in- 
teractions, and it is hard to pin down which effects are truly 
due to glassiness. Note that the somewhat simpler Bose glass 
of the Bose Hubbard model [18| (see [19] for a possible re- 
alization in optical cavities), while in the quantum regime, 
occurs because of a short-range random potential, and does 
not generically exhibit some of the hallmark phenomena of 
frustrated glasses, such as many metastable states, aging, or 
replica-symmetry breaking. 

It would clearly be desirable to have a tunable realization 



of genuinely frustrated (quantum) glasses in the laboratory. 
Recent work on ultracold atoms in optical cavities ll9l 4TTl[T9l 
suggests that it may be possible to create spin- and charge 
glasses in these systems, which arise because of frustrated 
couplings of the atomic "qubits" to the dynamical potential 
of multiple cavity modes. It is appealing to these systems 
that the photons escaping the cavity can be used for in-situ 
detection of the atom dynamics ("cavity glass microscope"), 
and that the interaction mediated by cavity photons is long- 
ranged. The latter makes the theoretical glass models more 
tractable and should allow for a realistic comparison of exper- 
iment and theory. 

The phases of matter achievable with cavity quantum 
electrodynamics (QED) systems settle into non-equilibrium 
steady states typically balancing a laser drive with dissipation 
channels such as cavity photon loss and atomic spontaneous 
emission. The notion of temperature is, a priori, not well de- 
fined. A line of recent research on the self-organization tran- 
sition of bosonic atoms in a single-mode optical cavity (ex- 
perimentally realized with a thermal gas of Cesium [20 1 and 
with a Bose-Einstein condensate of Rubidium lETl l22ll ). has 
established the basic properties of the non-equilibrium phase 
transition into the self-organized, superradiant phase [23- 32 J. 
In particular it was shown that, upon approximating the atom 
dynamics by a single collective spin of length N/2 and taking 
the atom number N large, the dynamics can be described by 
classical equations of motion [ 25 , 29 1 , and that the phase tran- 
sition becomes thermal due to the drive and dissipation [30|. 

In this paper, we underpin our previous proposal [101, an d 
show that quenched disorder from multiple cavity modes, 
leads to qualitatively different non-equilibrium steady states 
with quantum glassy properties. We develop a comprehen- 
sive non-equilibrium theory for many-body multimode cav- 
ity QED with quenched disorder and Markovian dissipation. 



We pay special attention to the quantum optical specifics of 
the pumped realization of effective spin model [ 33 1, the laser 
drive and the finite photon lifetimes of cavity QED. Using a 
field theoretic Keldysh formalism adapted to quantum optics, 
we compute several observables of the glass and superradi- 
ant phases, which are accessible in experiments with current 
technology. Our key results are summarized in the following 
section. 

The remainder of the paper is then organized as follows. 



In Sec. Ill we discuss the multimode open Dicke model in 
the simultaneous presence of quenched disorder and Marko- 
vian dissipation. Disorder and dissipative baths are contrasted 
more rigorously in App. [Bj We switch to a unified description 



of both these aspects in Sec. IV in the framework of a Keldysh 



path integral formulation, and specify the formal solution of 
the problem in the thermodynamic limit in terms of the par- 
tition sum, which allows to extract all atomic and photonic 
correlation and response functions of interest. This solution is 
evaluated in Sec.|V] with some details in App. [D] This com- 
prises the calculation of the phase diagram in the presence 
of cavity loss, as well as the discussion of correlation and 
response functions for both atomic and photonic degrees of 
freedom, allowing us to uniquely characterize the simultane- 
ous spin and photon glass phase from the theoretical perspec- 
tive. We then discuss the consequences of these theoretical 
findings to concrete experimental observables in cavity QED 
experiments in Sec. |VE| The combination of correlation and 
response measurements allows for a complete characteriza- 
tion of the phase diagram and in particular of the glass phase. 
The relation between Keldysh path integral and quantum 
optics observables is elaborated on further in App. IC] 



II. KEY RESULTS 




FIG. 1. (Color online) Non-equilibrium steady state phase diagram 
of the open multimode Dicke model, as a function of averaged atom- 
photon coupling J (y-axis) and disorder variance K (x-axis) and for 
parameters oj = 1 (cavity detuning) and to- = 0.5 (effective atom de- 
tuning) for different photon decay rates k. QG is the quantum spin 
and photon glass, SR the superradiant phase. The T = equilib- 
rium phase diagram of Ref. [ 10 1 is recovered as k — > 0. The SR-QG 
transition is not affected by k. 
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A. Non-equilibrium steady state phase diagram 

The shape of the phase diagram for the steady state pre- 
dicted in MOV , with the presence of a normal, a superradiant, 
and a glass phase is robust in the presence of Markovian dis- 
sipation, cf. Fig. U\ As to the phase diagram, the open nature 
of the problem only leads to quantitative modifications. In 
particular, the characteristic feature of a glass representing a 
critical phase of matter persists. The presence of photon de- 
cay overdamps the spin spectrum and changes the universality 
class of the glass phase, which we now discuss. 



B. Dissipative spectral properties and universality class 



FIG. 2. (Color online) Illustration of the dissipative spectral prop- 
erties and universality class. As a function of probe frequency u> (y- 
axis) and the disorder variance K (x-axis), we illustrate the different 
regimes in the phase diagram. In the normal phase, for frequencies 
a> < a the system is represented by a dissipative Ising model, de- 
scribed by Eq. (T21, while for frequencies co > a, uj c it is described by 
non-universal behavior of a disordered spin fluid. In the glass phase 
(K > K c ), there exist two qualitatively distinct frequency regimes, 
separated by the crossover scale co c , cf. Eq. {TJ. At the lowest fre- 
quencies, a> < oj c the system is described by the universality class of 
dissipative spin glasses. For cj > 0) c , we find that the system behaves 
quantitatively as an equilibrium spin glass. For a < co c and K < K c , 
there exists a dissipative crossover region (D-C in the figure), which 
is a precursor of the dissipative spin glass. It shows dissipative Ising 
behavior for smallest frequencies and resembles the dissipative glass 
for frequencies uj c > a> > a. 



Within the glass phase, we identify a crossover scale co c ~ k 
proportional to the cavity decay rate k, above which the spec- 
tral properties of a zero temperature quantum spin glass are 
reproduced. Although the finite cavity decay k introduces a 
finite scale "above the quantum critical point of the closed, 
equilibrium system", k acts very differently from a finite tem- 
perature. In particular, below co c , the spectral properties are 
modified due to the breaking of time reversal symmetry by the 



Markovian bath, while remaining critical. Due to the low fre- 
quency modification, the quantum spin glass in optical cav- 
ities formally belongs to the dynamical universality class of 
dissipative quantum glasses, such as glasses coupled to equi- 
librium ohmic baths [34-36] or metallic spin glasses 47711571/ . 
Spectral properties - The role of the crossover scale be- 
tween equilibrium and dissipative spin glass is further illus- 
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FIG. 3. (Color online) Dissipative spectral properties and univer- 
sality class of the single-atom spectral density ^H(a>) (response sig- 
nal of RF spectroscopy) in the quantum glass phase for parameters 
K = 0.01,7 = 0.1, a) z = %K = 0.1, 0>o = 0.7. For frequencies a> < io c 
below the crossover scale, the spectral density is overdamped and 
proportional to ^Joj. For intermediate frequencies a> > o) c , 3K is lin- 
ear in the frequency, as for the non-dissipative case 1101 . which is 
recovered in the limit k — > 0. 



trated in Fig. [2] It is given by 
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The resulting modifications below this scale, compared to a 
more conventional equilibrium glass are due to the Marko- 
vian bath, introducing damping. In the normal and superra- 
diant phases, this allows for the following form of the fre- 
quency resolved linearized Langevin equation for the atomic 
Ising variables, 



\ [to 2 + iyto + a 2 ) x(co) - t;(to), 



(2) 



modelling the atoms as an effective damped harmonic oscilla- 
tor, with finite life-time r = - < oo and a the effective oscil- 
lator frequency, with the physical meaning of the gap of the 
atomic excitations in our case. The noise has zero mean and 

<£(?')£«> = | reffc5(f-0. 

At the glass transition, Z and a scale to zero simultane- 
ously and the frequency dependence becomes gapless and 
non-analytic. In the entire glass phase, the effective atomic 
low frequency dynamics is then governed by the form 



i Vw 2 + y\co\ x(oj) = f(<y), 



(3) 



which obviously cannot be viewed as a simple damped oscil- 
lator any more. The broken time reversal symmetry manifests 
itself in y, y > 0, thus modifying the scaling for to — > 0. The 
crossover between these different regimes is clearly visible in 
Fig.|3] 

Universality class - The qualitative modification of the 
low-frequency dynamics below the crossover scale to c implies 
a modification of the equilibrium quantum spin glass univer- 
sality class. The open system Dicke superradiance phase tran- 
sition, where the Z2 symmetry of the Dicke model is broken 
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FIG. 4. (Color online) Thermalization into quantum-critical regime 
of the atomic (red, dashed line) and photonic (blue lines) distribution 
functions F(cS) when approaching the glass transition at a critical 
disorder variance K c for a> = 1.3, io z = 0.5, k = 0.01, K c = 0.01, 
J = 0.1 and varying parameter S = K c - K. For larger values of 
J, i.e. larger distance from the glass transition, the low frequency 
effective temperature (LET) 2T cif = lim^o oiF(u) of the photons is 
much lower than the LET of the atoms and the frequency interval 
for which atoms and photons are not equilibrated is larger. When 
the glass transition is approached, atoms and photons attain the same 
LET. 



spontaneously due to a finite photon condensate, is enclosed 
by a finite parameter regime in which the dynamics is purely 
dissipative, or overdamped (see e.g. [30 1). Together with the 
generation of a low frequency effective temperature (LET), 
for this reason the single mode Dicke phase transition can be 
classified within the scheme of Hohenberg and Halperin ||38l 
in terms of the purely relaxational Model A, thereby sharing 
aspects of an equilibrium dynamical phase transition. This sit- 
uation is different for the open system glass transition: Here, 
irreversible dissipative and reversible coherent dynamics rival 
each other at the glass transition down to the lowest frequen- 
cies. In particular, the dissipative dynamics fades out faster 
than the coherent dynamics as witnessed by larger critical ex- 
ponents, and there is no regime in the vicinity of the criti- 
cal point where either dissipative or coherent dynamics would 
vanish completely. This behavior is demonstrated in the inset 
of Fig. [7] 

We note that, while these findings are unconventional from 
the viewpoint of equilibrium quantum glasses, they are not 
uniquely tied to the presence of the driven, Markovian non- 
equilibrium bath. In fact, such behavior is also present in the 
case of a system-bath setting in global thermodynamic equi- 
librium, where the presence of the bath variables modifies the 
spectral properties of the spins ifTTl l34l [35). Both physical 
contexts share in common the time reversal breaking of the 
subsystem obtained after elimination of the bath modes and 
may be seen to belong to the same universality class. 



C. Atom-photon thermalization into quantum-critical regime 

As in the driven open Dicke model, the statistical properties 
of atoms and photons are governed by effective temperatures 
at low frequencies. The effective temperature differs in gen- 
eral for the two subsystems. Approaching the glass transition, 
these effective temperatures are found to merge. The finite 
cavity decay enables this mechanism but k does not directly 
play the role of effective temperature. This mechanism pushes 
the hybrid system of atoms and photons in the glass phase 
into a quantum-critical regime described by a global effective 
temperature for a range of frequencies. This quantum critical 
regime retains signatures of the underlying quantum critical 
point. 

The Markovian bath not only affects the spectral properties, 
but also governs the statistical properties of the system. The 
main statistical effect is the generation of a LET for the atomic 
degrees of freedom, for which we find 
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throughout the entire phase diagram, and taking the same 
value as in the single-mode case (in the absence of sponta- 
neous emission for the atoms). This thermalization of the 
atoms happens despite the microscopic driven-dissipative na- 
ture of the dynamics, and has been observed in a variety of 
driven open systems theoretically [28 39-45 1 and experimen- 
tally [46|. Below this scale, the occupation properties are gov- 
erned by an effective classical thermal distribution 2T e s/a), 
while above it the physics is dominated by non-equilibrium 
effects |30|. For cavity decay k <k ojq, the crossover scale 
obeys co c <K r dl . As a consequence, in an extended regime 
of frequencies between a> c and T e g, the correlations describe 
a finite temperature equilibrium spin glass. 

In the single-mode open Dicke model, the photon degrees 
of freedom are also governed by an effective temperature, 
which however differs from the one for atoms [30], indicat- 
ing the absence of equilibration between atoms and photons 
even at low frequencies. The increase of the disorder variance 
leads to an adjustment of these two effective temperatures, cf. 
Fig. HI At the glass transition, and within the entire glass 
phase, the thermalization of the subsystems is complete, with 
common effective temperature given in Eq. (HJ). This effect 
can be understood qualitatively as a consequence of the dis- 
order induced long ranged interactions, cf. Sec. ( |45] l. These 
allow to redistribute energy and enable equilibration. 

We emphasize that the notion of thermalization here refers 
to the expression of a 1 ju> divergence for the system's distri- 
bution function, as well as the adjustment of the coefficients 
for atoms and photons. This provides an understanding for 
distinct scaling properties of correlations (where the distribu- 
tion function enters) vs. responses (which do not depend on 
the statistical distribution), which can be addressed separately 
in different experiments (see below). Crucially, this notion of 
"thermalization" does not mean that the characteristics fea- 
tures of the glass state are washed out or overwritten. 
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FIG. 5. (Color online) Emergent photon glass phase with alge- 
braically decaying photon correlation function g (2 \r) at long times, 
for parameters cj = 1, k = 0.4, a>- = 6, J = 0.4, K = 0.16. The time- 
scale for which algebraic decay sets in is determined by the inverse 
crossover frequency co c , given by Eq. (75). For comparison, we have 
also plotted the envelope of the exponential decay of the correlation 
function in the normal and superradiant phase. The short time be- 
havior of the correlation function is non-universal and not shown in 
the figure, however, g <2) (0) = 3 due to the effective thermal distribu- 
tion for low frequencies. The parameter To = 0(:jr) was determined 
numerically. 



D. Emergent photon glass phase 

The strong light-matter coupling results in a complete im- 
print of the glass features of the atomic degrees of freedom 
onto the photons in the cavity. We refer to the resulting state 
of light as a photon glass highlighting the connection ofmul- 
timode cavity QED to random losing media l\47\ \48V . 

The photon glass is characterized by a photonic Edwards- 
Anderson order parameter signaling infinitely long memory 
in certain temporal two-point correlation function. This im- 
plies that a macroscopic number of photons is permanently 
present in the cavity (extensive scaling with the system size), 
which are however not occupying a single mode coherently, 
but rather a continuum of modes. The presence of a contin- 
uum of modes at low freqeuency is underpinned by the slow 
algebraic decay of the system's correlation functions as shown 
for the photon correlation function in Fig. [5] This is a conse- 
quence of the disorder-induced degeneracies. g (2) (r) is acces- 
sible by detecting the photons that escape the cavity. 



E. Cavity glass microscope 

The cavity set-up of Fig. \fi\should allow for unprecedented 
access to the strongly coupled light-matter phase with disor- 
der. Adapting the input-output formalism of quantum optics 
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FIG. 6. (Color online) Cavity glass microscope set-up: Atoms are 
placed in a multimode cavity subject to a transversal laser drive with 
pump frequency co p . The atoms are fixed at random positions by 
an external speckle trapping potential over regions inside the cavity, 
wherein mode functions g(k,, x/) randomly change sign as a function 
of the atomic positions, in order to provide frustration, as well as vary 
in magnitude. The more cavity modes, the better, and in particular 
the regime where the ratio of the number of cavity modes (M) over 
the number of atoms (N), a = M/N is kept sizable is a promising 
regime for glassy behavior |9. 15]. Photons leaking from the cavity 
with rate k give rise to additional dissipative dynamics and allow for 
output detection measurements. 



| \50V to the Keldysh path integral, we provide a compre- 
hensive experimental characterization of the various phases 
in terms of the cavity output spectrum and the photon corre- 
lations g™(T) in the real time domain. 

This continues and completes a program started in [30| of 
setting up a translation table between the language and ob- 
servables of quantum optics, and the Keldysh path integral 
approach. The frequency and time resolved correlations can 
be determined via fluorescence spectroscopy, cf. Fig. [71 and 
the measurement of g (2 \r) follows time-resolved detection of 
cavity output, cf. Fig. 15] The fluorescence spectrum shows 
a characteristic -4= divergence for small frequencies to < co c . 
This indicates a macroscopic but incoherent occupation of the 
cavity as anticipated above: The glass state is not character- 
ized by a single-particle order parameter where a single quan- 
tum state is macroscopically occupied, and which would re- 
sult in (temporal) long range order such as a superradiant con- 
densate. Rather it is characterized by a strong and infrared 
divergent occupation of a continuum of modes, giving rise to 
temporal quasi-long range order. This phenomenology is rem- 
iniscent of a Kosterlitz-Thouless critical phase realized e.g. in 
low temperature weakly interacting Bose gases, with the dif- 
ference that spatial correlations are replaced by temporal cor- 
relations. 

Finally, the combined measurement of response and corre- 
lations enables the quantitative extraction of the effective tem- 
perature. 
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FIG. 7. (Color online) Cavity glass microscope output of a typ- 
ical fluorescence spectrum S(a>) (not normalized), decomposed in 
coherent S c and incoherent part 5 inc for the three distinct phases in 
the multimode Dicke model. The parameters J, K are varied, while 
a> = 1,k = 0.1,cl>- = 0.5 are kept fixed for each panel. 
Normal phase, (J, K) = (0.13, 0.008). Central and outer doublets are 
visible but broadened by the disorder, only the incoherent contribu- 
tion is non-zero. 

Superradiant phase, (J, K) = (0.4, 0.008). The central doublets have 
merged due to the presence of a critical mode at CO = 0. At zero 
frequency there is a coherent ^-contribution indicated by the arrow 
(dashed). 

Glass phase, (J,K) = (0.13,0.017). There is a characteristic -4= di- 
vergence for small oj < uj c due to the non-classical critical modes at 
zero frequency. The peak at a> = is incoherent and can therefore 
easily be discriminated from the coherent peak in the middle panel. 
Scaling of correlations. The inset in the upper panel shows the be- 
havior of the peak distance of 5 (to) in the normal phase when ap- 
proaching the glass phase. The two peaks approach each other and 
merge at the glass transition. The dista nce follows the dominant co- 
herent exponent a 6 oc £5 , cf. Sec. II B 



III. MULTIMODE OPEN DICKE MODEL 

In this section, we explain the model for fixed atoms in 
an open multimode cavity subject to a transversal laser drive 
shown in Fig. [6] We first write down the explicitly time- 
dependent Hamiltonian operator for a level scheme involving 
two Raman transitions proposed by Dimer et al. [33). We 
then transform this Hamiltonian to a frame rotating with the 
pump frequency. This eliminates the explicit time dependen- 
cies in the Hamiltonian at the expense of violating detailed 
balance between the system and the electromagnetic bath sur- 



rounding the cavity. The bath becomes effectively Markovian 
in accordance with standard approximations of quantum op- 
tics. Finally, we eliminate the excited state dynamics to arrive 
at a multimode Dicke model with tunable couplings and fre- 
quencies. 



A. Hamiltonian operator 

The unitary time evolution of the atom-cavity system with 
the level scheme of Fig.JHlfoilows the Hamiltonian 
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which we now explain one-by-one. The cavity contains M 
photon modes with frequencies v,- 



|o> 



Z/cav - Zj=l ViQ-Qi, 



(6) 



which we later take to be in a relatively narrow frequency 
range v, a vo such that the modes couple with comparable 
strengths to the detuned internal transition shown in Fig. [8] 
The atom dynamics with frequencies given relative to the 
lower ground state |0) is 



FIG. 8. (Color online) Internal level scheme to generate tunable 
Dicke couplings between the ground state levels |1), |0> and the cav- 
ity. Adapted from Dimer et al. 1331 . 



with a correspondence of the effective spin operators in 



Eq. ( 10 1 to the internal atomic levels 



o\ = |1*><1<I - |0*>«M , cr£ = |1,><0<| + 10/XWI . (11) 



// a t - y. 0)rVt){rt\ + a> s \sc){s[\ + a»i|l^)(lf|. (7) The couplings and frequencies are tunable: 



The interaction between the atoms and cavity modes 

N M 

//i„, = Yj Z (&-(k*,*«)l»v>«M + g,Qn,*i)\s£>(Utyh + H.c. 

(8) 

involves a set of cavity mode functions g(k,, Xf) which depend 
on the wave vector of the cavity mode k, and the position of 
the atom Xf , The pump term 



0)i = V; -(w p>s - Wj) + 

g r (ki,x ( )Q. r 
8it ~ 2A r ' 



~gf(h) 



2(wi — u>\), 



(12) 



2 
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does not involve photon operators and induces coherent tran- 
sitions between the excited and ground states as per Fig. [8] co p 
is the (optical) frequency of the pump laser. We assume the 
atoms to be homogeneously pumped from the side so that the 
mode function of the pump lasers are approximately constant 

We now transform Eqs. (|5]l-((9]l to a frame rotating with the 
(optical) frequency of the pump laser, mainly to eliminate the 
explicit time dependence from the pump term Eq. (|9]) [ 33 1 . 
The unitary transformation operator is U(t) = exp(— ifltf) 
with H = [a>p, s - a>\) Z^i a\a t + J$ ml { (w p>r + oj\) \r c ){r c \ + 

<a>p,s\s{){S{\ + w' 1 |l^)(lf[J , with u)\ a frequency close to o>i 
satisfying a> PtS - u> pj - = 20*^ ll33l . We then eliminate the ex- 
cited states in the limit of large detuning A to finally obtain 
the multimode Dicke model 

h = zf =l u$ai + t Efei °\ + Zu f^/ v ( S J + Si) , (io) 



.In particular the effective spin-photon cou- 
pling gic can now be tuned sufficiently strong to reach superra- 
diant regimes by changing the amplitude of the pump Q r . The 
effective cavity frequencies receive an additional shift from a 
mode mixing term a,flj with space averaged cavity couplings 
~ g 2 /A r from which we only keep the mode-diagonal contri- 
bution (for running wave cavity mode functions ~ e' k,Xf this 
is exact; we do not expect qualitative changes to our results 
from this approximation). 

The multimode Dicke model with internal atomic 
levels obeys the same Ising-type Z2 symmetry, 
(a l , cr*) — > (—a. , -erf), familiar from the single-mode 
Dicke model fl5"flj56ll . Therefore, there exists a critical 
coupling strength J c , such that the ground state of the system 
spontaneously breaks the Z2 symmetry as soon as the average 
coupling strength 
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exceeds the critical value, J > J c , The phase transition from 
the symmetric to the symmetry broken, superradiant (SR) 
phase has been well analyzed for the single-mode Dicke 
model and the essential findings, such as the universal behav- 
ior for zero and finite temperature transitions 1 55 56 1 or in the 
presence of dissipation [27 , 28 30], remain valid also for the 
multimode case. The superradiant phase is determined by the 



presence of a photon condensate, i.e. the emergence of a co- 
herent intra-cavity field [22 57 1, which is described by a finite 
expectation value of a photon creation operator (a c ) + 0. The 
superradiant condensate a c - Yji a< i a > with Z,-|a^| 2 - 1. 
is a superposition of many cavity modes at, and its explicit 
structure depends on the realization of the couplings {gu}. 



B. Markovian dissipation 

In a cavity QED experiment of the type described in Fig.|6j 
the atoms and photons governed by the Hamiltonian ( | 1 Op 
are additionally coupled to the electromagnetic field outside 
the cavity. This leads to the additional processes of sponta- 
neously emitted photons into the environment and to cavity 
photon loss through imperfect mirrors, accurately captured by 
a Markovian master equation |58 , 59 1) of the form 



d t p = -i[H,p]+£(p) = £ tat (p), 



(14) 



where p is the density matrix of the atom-photon system, H is 
the Hamiltonian (jTOJ and X is a Liouville operator in Lindblad 
form 



-C(f) = Ha K a (2L a pLl - {LlL a ,p}J . 



(15) 



Here, {■ ,} represents the anti -commutator and the L a are 
Lindblad or quantum jump operators. The photon dissipation 
is described by the Liouvillian 



-C P h(p) = Zfei Ki (2a,-pa,- - {a\ai,p}) , 



(16) 



where k, is the loss rate of a cavity photon from mode (/)■ 
Eq. ([16) describes a Markovian loss process that, while being 
a standard approximation in quantum optics, violates detailed 
balance between the system and the bath. Formally, it can be 
derived by starting with a cavity-bath setup in which both are 
at equilibrium with each other, and performing the transfor- 



mation into the rotating frame outlined above Eq. ( 10 1 also on 



the system-bath couplings (see App. B 3 1. 

In this work, we consider Ki < U)i, u> : but of the same or- 
der of magnitude. In contrast, the atomic dissipative dynam- 
ics are considered to happen by far on the largest time scale, 
which can be achieved in typical cavity experiments ll22l 1571 . 
In a recent open system realization of the single-mode Dicke 
model |2"2"ll57l . spontaneous individual atom-light scattering 
is suppressed by five orders of magnitude compared to the rel- 
evant system time-scales, such that atomic dephasing effec- 
tively plays no role [22J. As a result, only global atomic loss 
is influencing the dynamics, which, however, can be compen- 
sated experimentally by steadily increasing the pump inten- 
sity or chirping the pump-cavity detuning [22] . We therefore 
do not consider atomic spontaneous emission in this paper. 



is sufficiently large. The specific values of the couplings gu 
in Eq. ( fTO) are fluctuating as a function of the atom (/) and 
photon (/) numbers and depend on the cavity geometry and 
realization of the random trapping potential (Fig. 15). After 
integrating out the photonic degrees of freedom in Eq. ( [TO) , 
we obtain the effective atomic Hamiltonian 



H M 






2w,,„=l JlmO~iO-„ 



(18) 



where we introduced the effective atom-atom couplings J\ m = 
2^i ^|^-, and at this point neglected the frequency depen- 
dence in the atom-atom coupling term in Eq. ( fT8) . This is ap- 
propriate for a>, « a> and wo large compared to other energy 
scales (in particular, |(w, - a>o)/a>o\ <K 1). In order to solve 
the effective Hamiltonian ( p"8) , it is sufficient to know the dis- 
tribution of the couplings // m , which itself is a sum over M 
random variables. For a large number of modes (M — > oo), 
this distribution becomes Gaussian, according to the central 
limit theorem, with expectation value J and variance K, as 
defined in Eqs. ( |13) , ( fP7) , respectively. 

The variables //,„ can be seen as spatially fluctuating but 
temporally static variables, connecting all atoms with each 
other. This may be seen as a coupling to a bath with random 
variables J\ m , which vary on time scales tq much larger than 
the typical time scales of the system t s only. The dynamics 
of the bath is then frozen on time scales of the system, and the 
bath is denoted as quasi-static or quenched [14|. This type of 
bath is in a regime opposite to a Markovian bath, where the 
dynamics of the bath happens on much faster time scales t m 
than for the system, tm ^ ts [58, 59 1. We have summarized 
basic properties of these baths in App.[B| 



IV. KELDYSH PATH INTEGRAL APPROACH 

In this section, we introduce the Keldysh formalism [30 
1501 I5T1 and derive the set of self-consistency equations for 
the atoms and photons from which all our results can be ex- 
tracted. We first formulate the open multi-mode Dicke model 
Eqs. ( 10|16 l as an equivalent Keldysh action that includes the 
non-unitary time evolution induced by cavity decay. In the 
Keldysh approach, one additionally benefits from the fact that 
the partition function 



Z = Tr(pO0) = l 



(19) 



is normalized to unity, independent of the specific realiza- 
tion of disorder, and we perform the disorder average directly 
on the partition function. We then integrate out the photons 
(carefully keeping track of their correlations, as explained be- 
low) and derive a set of saddle-point equations for frequency- 
dependent correlation functions which can be solved. 



C. Quenched / quasi-static disorder 



A. Multi-mode Dicke action 



The glassy physics addressed in this paper arises when the 
spatial variation of cavity mode couplings 

N 2 



Y - 1 V" (V M Siigim Y _ j2 
"■ — N t-il,m=\ \lji=l 4 ) J 



(17) 



To describe the photon dynamics, one starts from an ac- 
tion for the coupled system of cavity photons and a Marko- 
vian bath. Then the bath variables are integrated out in Born- 
Markov and rotating wave approximations. The resulting 
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Markovian dissipative action for the photonic degrees of free- 
dom on the (±)-contour reads 

dt(a* + (id, - Mj)a j+ - (a*_(id, - iOj)a.j_ 

co 



- iK[2a. + a*_ 



- (a j+ a j+ + a). 



ijjl). 



(20) 



Here, the creation and annihilation operators have been re- 
placed by time-dependent complex fields. The structure of 
the master equation (|T4bis clearly reflected in the action on 
the (±)-contour in Eq. (|20|i. The first line corresponds to the 
Hamiltonian part of the dynamics, with a relative minus sign 
between (+) and (-) contour stemming from the commutator. 
The second line displays the characteristic form of the dissi- 
pative part in Lindblad form. 

For practical calculations, it is more convenient to 
switch from a (±)-representation of the path integral 
to the so-called Keldysh or RAK representation. In 
the latter, the fields on the (±)-contour are transformed 
to "classical" a jc = (o + +a_)/y2 and "quantum" fields 

a . — (a . + — a ._)/ v2, where the labeling of these fields in- 
dicates that a . can acquire a finite expectation value, while 
the expectation value of a . is always zero. After a Fourier 

transformation to frequency space, a t (co) - J dt a t (t) e~'°", 
the photonic action in Keldysh representation is obtained as 



i 



,W*c 



«},P 



D R (co) 
D A (co) Df(co) 



a . 

M 



(21) 



where we used the abbreviation J. = 2 %i J ff • The integral 
kernel of Eq. pT| ) is the inverse Green's function in Keldysh 
space with the inverse retarded/advanced Green's function 



DfV) = [Gf A r» 



CO + IK; 



(22) 



and the Keldysh component of the inverse Green's function 

Df(co) = 2k,-. (23) 

j j 

From now on, we will focus on the case where the variation 
in the photon parameters is much smaller than all other en- 
ergy scales of this problem and consider only a single pho- 



ton frequency coq and photon loss rate k, i.e. 



<k k and 



\(Oq - U)j\ 



<§c cjq for all photon modes (j). As a result all pho- 



ton Green's functions are identical with Kj - k and ojj 
The Green's function in Keldysh space takes the form 



0{oS) 



G K (u) G R (to) 
G A (u>) 




D A (u>) 



D R {cS) 
D K (to) 



loo. 



(24) 



where we already identified retarded/advanced Green's func- 
tion G R {a>) in Eq. ( f2"2"| ). The Keldysh component of the 
Green's function is obtained by performing the inversion (|24j» 

as 



G K (co) = -G r (cj)D k (co)G a (oj). 



(25) 



The retarded Green's function encodes the response of the sys- 
tem to external perturbations and its anti-hermitian part is pro- 
portional to the spectral density 



m(io) = i(G R (co) - G A (u>j) , 



(26) 



since G A (a>) - \G r {(jS)\ . The retarded Green's function 

G r (cl>) and the Keldysh Green's function G K (u>) constitute the 
basic players in a non-equilibrium path integral description, 
determining the system's response and correlations. For a 
more detailed discussion of a Keldysh path integral descrip- 
tion of cavity photons, we refer the reader to 1 30 1 . 

The atomic sector of the Dicke Hamiltonian ( fLO] ) can be 
mapped to an action in terms of real fields (pi, as long as the 
physically relevant dynamics happens on frequencies below 
co : |62|. The <f>[ obey the non-linear constraint 



6(tf(t) - 1) = fDAiity- 



1,(0(^(0-1) 



(27) 



where Ai(t) are Lagrange multipliers, in order to represent 
Ising spin variables (see Ref. 30 for further explanation). As 



a result, we can apply the following mapping to Eq. ( 10 1 



<Pi(t), 



1, 



(28) 
(29) 



<r\(t) 

LL/ 7 

On the (±)-contour, we then obtain 

S«=ih t (dt<Pi + f-(d t <!>i-) 2 , (30) 

subject to the non-linear constraint 

s~ = £.fcM&-i)-M£-i). (31) 

The atom-photon coupling reads 

S„ UP = X„ f (</>,+ (a! + + a i+ ) - fr. (al + «._)) . (32) 

Transforming to the RAK basis and frequency space, the 
atomic propagator becomes 

S m = — f {fa,fa)DJfii)l* c * )+— f V (33) 



with 



DJfa) = 



A cJ - (oj - ir/) 2 



v,/ 



(o» + irj) 



I 



q,l 



(34) 



Here, r\ — > + plays the role of a regulator that ensures causal- 
ity for the retarded/advanced Green's functions. For the atom- 
photon coupling in the RAK basis, we get 



S„ up = I Y \(<PcJ,(f>qj)o- X 



q,l ) 



. / * * \ x 

+ ( a C.l> a ClP Cr 



4>C,l 

<t>q,l 



(35) 



For the atomic fields, it is useful in the following to introduce 

0c/(w) 



the Keldysh vector <3>/(<y) = 



(/>q,i(<o) 



, which will simplify 



the notation in the following. 

The Keldysh action for the open multimode Dicke model is 
then obtained as the sum of Eqs. ( 21|33 35 1 



S[{a f ,a ,<t>,A}] 



" ph + & al + " coup • 



(36) 



Atoms 


Photons 








G«ftO = Q K (t,f) = -i («(?),<(?')}) 

&»(*■ O = Q R (t, O = -i ®(f - O ([cr* (t), of (C)]) 

Q qq (t,f) = 


G CC (U') = G*(f,0 = -i (K(0,al(f)}) 

G cg (f, = G R (f, t) = -i &(t - f) ([a,„(t), al(t' )]) 

G 9C (f,f) = G*(.t,f) = (G*(t,fjf 

G w (f,O = 


<M0 = V2 <of (r)> 

<A,M = 



TABLE I. Translation table for the atomic order parameter and Green's functions, from now on labeled with Q, and the intra-cavity photon 
Green's function, labeled with G 



B. Calculation procedure 



We now explain how we solve the Keldysh field theory de- 



value. The disorder averaged partition function can be 
expressed as 



scribed by Eq. (36 1. The calculation proceeds in three steps: 



1. Integration of the photon modes: This step can be per- 
formed exactly via Gaussian integration, since the action |36]l 
is quadratic in the photon fields. Note that this does not mean 
that we discard the photon dynamics from our analysis. To 
also keep track of the photonic observables, we modify the 
bare inverse photon propagator, Eq. ( [21) , by adding (two- 
particle) source fields yU according to 



Z = JD({®,A,J})e i{S « +s ^ +s ^\ 
with the disorder "action" 



(43) 



D ph (w) -> D^to) + //(w), with fi 



p? c p qq 



(37) 



The photon Green's functions are then obtained via functional 
variation of the partition function with respect to the source 
fields 



G*' A > K (a>) 



Spfl'Mcc^) |^ =() ■ 



(38) 



The resulting action is a sum of the bare atomic part ( |33) and 
an effective atom-atom interaction 



Sdis - 2 2jl,mJ',m' (/lm Jim) ^bnl'm' {Jl'tri Jl'm')> (44) 

describing a temporally frozen bath with variables J[ m . Per- 
forming the disorder average, i.e. integrating out the variables 
Jim in the action ( |43"| > replaces the parameters 7; m — > J/N in 
(|39j) by their mean value. Furthermore, the variance K intro- 
duces a quartic interaction term for the atomic Ising variables 
which is long-range in space 

S,, = f JL, Zi, m (a»/A(D m ) (oS) (*,AO m ) (o/), (45) 

with the shortcut (0;A<J> m )(w) = ®J(a))A(a))® m (ci)), 

3. Collective variables: Atomic order parameter and 
Green 's function: To decouple the spatially non-local terms in 
([39]) and (J45J, we introduce the Hubbard-Stratonovich fields 
i/r„ and Q aa i with a, a' = c, q, which represent the atomic or- 
der parameter 



tffaiaj) = j, £/<0ff,/M> 



S« = - £ Zi,m Ji m ®J(-«))A(a>)® m (a)), (39) and avera 8 e atomic Green ' s function 



with atom-atom coupling constants //,„ defined in ( fl8] i and the 
frequency dependent coupling 



QaAoiM) = i Z/<<^,/M0„ V (o/)>. 



(46) 



(47) 



A(oj) = \cr x (G (o>) + Gj(-w)) o- x , 



(40) 



which is the bare photon Green's function Go after sym- 
metrization respecting the real nature of the Ising fields <£/. 
We note that the information of the photonic coupling to the 
Markovian bath is encoded in A(w). 

2. Disorder average: The coupling parameters J\ m are con- 
sidered to be Gaussian distributed and the corresponding dis- 
tribution function is determined by the expectation value and 
covariance of the parameters J/,„ 



Now, the action is quadratic in the original atomic fields <pt, 
and so these can be integrated out. The resulting action has 
a global prefactor N and we will perform a saddle-point ap- 
proximation which becomes exact in the thermodynamic limit 
and upon neglecting fluctuations of the Lagrange multiplier. 
We replace the fluctuating Lagrange multipliers Ai(t) by their 
saddle-point value Ai(t) = A. In the steady state, the atomic ob- 
servables become time-translational invariant which restricts 
the frequency dependence of the fields to 

ifr a (.to) = 2n5(oj)i/j a , (48) 

Q aa >(u>, o/) - 2tt5(u + co')Q aa ,(a>). (49) 



6Ji m 6Ji>, 



Jim ~ N , 

K 



N 



(6u'8mm' + <>lm' $ ml') = %l 



Iml'i 



(41) 
(42) 



where the line denotes the disorder average and 
5 Jim - Jim - Jim represents the variation from the mean 



C. Saddle-point action and self-consistency equations 

The saddle-point action is given by the expression 
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S/N 



2,1„ 



^ T (-co) [jA(to) - J 2 A(co)G(co)A(co)] ¥(w) - |Tr [lnG(w)] + jJTTr [AgAg] (w), 



(50) 



with the "Green's function" 

G(co) = (d m (w) - 2KA(co)Q(co)A(co)j (5 1 ) 

and the field W 7 = (iff c ,iff q ). The matrices A,G and 2 in 
Eq. |50) possess Keldysh structure, i.e. they are frequency 
dependent triangular matrices with retarded, advanced and 
Keldysh components. The matrix A contains the photon fre- 
quencies (l>o, the decay rate k, and also depends on the photon 
Lagrange multiplier p, so that all photon correlations can be 
extracted from Eq. d50b. 



1. Atomic sector 

In order to find a closed expression for the macroscopic 
fields {<t>, Q] and to determine the saddle-point value for the 
Lagrange multiplier A, we have to evaluate the saddle-point 
equations 



■>.v 



= 0, With X - Qaa',l//a, ^-a, a - C,q. 



(52) 



In stationary state, A q — iff q - Q qq — by causality and we set 
A c — A and i// c - iff for convenience. 

The saddle-point equation for A q expresses the constraint 



2 = fjQ K (co) = iQ K (t=0) 



2iZf=i<K) 2 >, 



(53) 



which has been reduced to a soft constraint, present on av- 
erage with respect to (I), compared to the original hard con- 
straint, (cr*) 2 - 1 for each (0 individually. 

In the superradiant phase and in the glass phase, the spin 
attain locally "frozen" configurations. The correlation time of 
the system becomes infinite, expressed via a non-zero value 
of the Edwards-Anderson order parameter 



<7e 



lim T _ 



W Z£l K(T)<7?(0)>. 



(54) 



As a consequence, the correlation function Q K (co) is the sum 
of a regular part, describing the short time correlations and a 
^-function at co — 0, caused by the infinite correlation time. 
We decompose the correlation function according to 



Q K {u) = 4inq^8((o) + Q K Jco), 



(55) 



with the Edwards-Anderson order parameter q EA , being de- 
fined in Eq. |54} and a regular contribution g£. In the lit- 
erature [34, 63 1, this decomposition is referred to as modi- 
fied fluctuation dissipation relation (FDR) as also discussed in 
Appendix IB] The saddle-point equations for atomic response 
function and the regular part of the Keldysh function are 

eV) = (^ - 4*(aV)) 2 eV))~' (56) 



and 



Gf,(w) 



4Ar|G R | 2 A s '(e' , A A +G R A*) 

l-4tf|e R A R | 2 



(57) 



Eqs. ( |53"j ), ( |56"l l, (J57j form a closed set of non-linear equations, 
describing the physics of the atomic subsystem in the thermo- 
dynamic limit, which will be discussed in Sec. IV] 



2. Photonic sector 



The photon response G R and correlation function G K are 
determined via functional derivatives of the partition function 
X, with respect to the source fields jj, as described in ( f3"7] l and 
(|38)). The saddle-point for the partition function is 



Z = e' b x Z 



r(0) 



(58) 



with the action S from Eq. |50} and the bare photon partition 
function Z ph . 

In the Dicke model, the photon occupation n-, is not a 
conserved quantity, such that anomalous expectation values 
(a 2 ) + will become important. This has to be taken into 
account by introducing a Nambu representation, where the 
photon Green's functions become 2x2 matrices, see Ap- 
pendix [A] Generalizing the source fields // to include normal 
and anomalous contributions, and evaluating the functional 
derivatives with respect to these fields, results in the inverse 
photon response function 



£?X2(W) - 



(59) 
(60) 



CO + IK ~ O> + £ R (w) £*(to) 

(S ft (-w))* -co - ik - oj + (E s (-w)) 

Here, the subscript 2x2 indicates Nambu representation and 

^) = (iV))*^(j||-i) 



(61) 



is the self-energy, resulting from the atom-photon interaction. 
The Keldysh component of the inverse Green's function is 

_ /2w+E*(a>) 2 K (fi>) \ 

D 2*2W = [_(xK ((0) y 2iK-(^(0S)JJ ^ 

with the self-energy 

v ' V v ') 4Re(e ff (w)A fl (w)) v ' 

In the Dicke model, the natural choice of representation for 
the photon degrees of freedom is the x-p basis, i.e. in terms 
of the real fields x — 4= (a* + a), p — -j=-.( a * ~ o), since 
the atom-photon interaction couples the photonic ^-operator 
to the atoms. In this basis, the self-energy gives only a contri- 
bution to the x-x component of the inverse Green's function, 
and the inverse response function reads 



<x>0 K — ICO 

-k + ico -COq 



(64) 



In the limit of vanishing disorder K — > 0, the self-energy ap- 



proaches the value £ R (w) = 



Joj, 



2{u> 2 -A 

the single mode Dicke model 11301155 



, reproducing the result for 
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V. RESULTS 

We now present our predictions from solving the atomic 
saddle-point equations Eqs. |53|, ([56]), ((57]) and then extract- 
ing the photonics correlations using Eqs. (|60] i-(|64|i, in the 
same order as in the Key Results Section |Il| In the subsec- 
tion Cavity Glass Microscope, we present signatures for stan- 
dard experimental observables of cavity QED by adapting the 
input-output formalism to the Keldysh path integral. 



A. Non-equilibrium steady state phase diagram 

The phases in the multimode Dicke model shown in Fig.[T] 
can be distinguished by the two order parameters, namely the 
Edwards-Anderson order parameters q EA in Eq. |54|, indicat- 
ing an infinite correlation time t and the ferromagnetic order 



parameter if/ defined (Eq. (48 1), indicating a global magneti- 
zation: 



Normal 


<7ea = 


«A = 


SR 


<7ea ^0 


tA*0 


QG 


<7ea * 


tA = 



In the normal phase, the Edwards-Anderson parameter q EA 
and the ferromagnetic order parameter ijj are both zero and 
Eq. <[53j implicitly determines the numerical value of the La- 
grange parameter A N . In contrast, in the superradiant phase 
if/ + 0, and the Lagrange parameter can be determined analyt- 
ically to be 



ASR - wJ+K 2 V + J)' 



(65) 



In the quantum glass phase the Lagrange parameter is pinned 
to 



'Iqg = 



V*. 



(66) 



The normal phase is characterized by a vanishing Edwards- 
Anderson order parameter, and the corresponding Lagrange 
multiplier A H is determined via the integral 



o = 2-/£e»| 



i-iN 



(67) 



The normal-SR phase border is located at the line for which 
A N = A SR , while the normal-QG transition happens at A N - A QG . 
In the same way, the transition between superradiant phase 
and quantum glass phase happens when \f/ vanishes for finite 
q EA + 0. This is the case for 



Ann <=> K = J 



(68) 



The phase diagram for the open system for different values of 
the photon dissipation k is shown in Fig. [T] As can be seen 
from this figure, the qualitative features of the zero tempera- 
ture phase diagram 1 10 1 are preserved in the presence of dis- 
sipation. However, with increasing k, the phase boundaries 
between normal and SR, QG phase are shifted to larger values 
of J, K respectively, while the QG-SR transition is still lo- 
cated at the values for which J 2 — K as for the zero tempera- 
ture equilibrium case. Finite dissipation neither favors the QG 



Analytic expressions 


SR to QG 


Normal to QG 


W 2 +K 2 

^fi 1 C T A 

" 8V^ 5 «r 2 


si 

s 2 
s 3 


i 

6 * 

log(i) 
i i2 

logW) 

1 * I' 
1 togOS) 1 



TABLE II. Atomic spectral response and scaling behavior approach- 
ing the glass transition in two different ways. The normal to glass 
transition shows logarithmic corrections compared to the SR to QG 
transition. The logarithmic scaling correction is a typical feature of 
the glass transition and has also been found for T = and finite tem- 
perature glass transitions in equilibrium |66 67 1. We see that the life- 
time of the excitations y 6 scales differently from the excitation energy 
as, which indicates a strong competition of the reversible quantum 
dynamics and the classical relaxational dynamics. Although the in- 
verse life-time scales faster to zero than the excitation energy, there 
is no point before the transition, where one of these quantities be- 
comes exactly zero as it was the case for the superradiance transition. 
The described behavior at the glass transition means that there is no 
purely relaxational low energy theory which is able to describe the 
dynamics close to the transition. It does not fall into the Halperin- 
Hohenberg classification of dissipative dynamical systems, but be- 
longs to the universality class of dissipative spin glasses 13444371 . 



nor the SR phase and as a result, the competition between dis- 
order and order is not influenced by the dissipative dynamics. 
The line at K — 0, i.e. zero disorder, describes the normal-SR 
transition for the single mode Dicke model, which is known to 

2 . „2 

be located at J c - |- a> z BUI IB31 l56l . This result is exactly 



reproduced within our approach. 



B. Dissipative spectral properties and universality class 

The atomic excitation spectrum and the influence of the 
system-bath coupling on the atomic dynamics are encoded 
in the retarded atomic Green's function, which is identical to 
the atomic linear susceptibility, Q R (u>) - %■ '(co). It describes 
the response of the atomic system to a weak perturbation as, 
for instance, the coupling to a weak coherent light field (see 



appendix C 1 1, and its imaginary part determines the atomic 
spectral response 



3{(u>) = -2lm(Q R (ojj), 



(69) 



which can be measured directly via radio-frequency spec- 
troscopy [64, 65 ) . The spectral response J[ for the normal and 
SR phase (the regular part for the latter) is shown in Fig. [9] 
In order to describe the low frequency behavior of the atomic 
spectrum, we decompose it into a regular and singular part, 
where the singular part captures the critical mode of the SR 
phase in terms of a ^-function at zero frequency, which is ab- 
sent in the normal phase. The regular part of the spectrum has 
the same structure for the normal and superradiant phase, and 
a derivative expansion of the inverse Green's function yields 
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the low frequency response function 

Q H (co)=Z s ((co + iy s ) 2 -a})~ 



(70) 



with the analytic expressions for the coefficients given in Ta- 
ble III] This is the Green's function of a damped harmonic os- 
cillator with characteristic frequency a> - as and damping y$, 
which is described by classical relaxational dynamics and cor- 
rectly determines the atomic spectral response for frequencies 
a> < \\ag - iy$\\ smaller than the pole. The index 6 in Eq. ( f70] i 
indicates that the parameters scale with the distance to the 
glass transition 



5 = K r -K, 



(71) 



which happens at K — K c (6 > in SR and normal phases). 
Transforming the response function to the time-domain, 



Q R (t) = Zse-*' cos(a s t) ., 



(72) 



which describes an excitation of the system with inverse life- 
time y$ = 1 /ts, energy as and quasi-particle residuum Zs. For 
frequencies a> < as, the spectral response is determined by the 
imaginary part of Eq. (|70b, yielding 



Mco) 



ZsJstii 



2k 



)VfiS 



CO. 



(73) 



This linear behavior is completely determined by parame- 
ters of the quenched and the Markovian bath and vanishes 
for k — * 0, resulting in a gap in the spectral weight for the 
zero temperature equilibrium case, as discussed in iflOl . For 
<jd > as, i.e. for a> larger than the gap, the atomic response 
function develops a non-analytic square root behavior 



J[{<x>) oc -\faJ~- 



as, 



(74) 



resulting from the quadratic form in Eq. ( |56) . The scaling 
of the excitation gap a 6 and the ratio -^ can be obtained 

directly from the atomic spectral response, as illustrated in 
Fig. [9] lower panel. At the glass transition, the gap vanishes, 
such that the square root behavior starts from u> - 0. 

In Table |H] we compare the scaling behavior of the atomic 
spectral response close to the normal-QG transition versus 
SR-QG transition lines. At the glass transition, Z s , js ar, d a s 
vanish, which for the latter two results in zero energy excita- 
tions with infinite life-time and therefore infinite correlation 
times. The vanishing of the residuum Zs indicates that the 
discrete poles of the system, representing quasi-particles with 
weight Z, transform into a continuum represented by a branch 



cut in the complex plane as illustrated in Fig. 10 As a conse- 
quence, a derivative expansion of the inverse propagator is no 
longer possible in the quantum glass phase. 

When approaching the glass phase, the frequency inter- 
val which is described by classical relaxational dynamics 
(i.e. [0, as]) shrinks and vanishes completely at the transition, 
where the system becomes quantum critical. The linear scal- 
ing of J[{<jS) in combination with the closing of the spectral 
gap is taken in thermal equilibrium as the defining property of 
a quantum glass. However, for a general non-equilibrium set- 
ting, the closing of the spectral gap is only a necessary but not 
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FIG. 9. (Color online) Regular part of the spectral density y{(oj) 
in the superradiant phase for parameters K = 0.05 and J = 0.4 and 
varying k and o>o- For small frequencies to < as the spectral den- 
sity is linear in a> and k and behaves as a square root for inter- 
mediate frequencies cj > a 6 . For the non-dissipative case {k — > 0), 
the spectral weight develops a gap at low frequencies, which is in- 
dicated for k = 10~ 3 (solid line). The lower panel depicts the low 
frequency behavior of ^H (red, dash-dotted line) for values k = 0.03 
and u) = 0.9. The green (full) and the black (dashed) line indicate 
the linear, square root behavior, respectively. Approaching the glass 
transition, a s scales to zero oc J 2 . 



a sufficient condition for the glass phase. The unique prop- 
erty of the glass transition in a non-equilibrium setting is the 
emergence of a critical continuum at zero frequency, which 
leads to the closing of the gap of the retarded Green's func- 
tion (distinct from the spectral gap). From the structure of the 
low frequency response function, Eq. ( fTO] ), we see that closing 
the spectral gap and a linear behavior of the spectral density 
is a non-trivial (and glass) signature only for a system where 
time-reversal symmetry is preserved, i.e. y — 0. On the other 
hand, the spectral gap closing is always present for broken 
time-reversal symmetry. 

Within the glass phase, it is again possible to separate 
two distinct frequency regimes delimited by a cross-over fre- 
quency 

which depends on all system and bath parameters. For 
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FIG. 10. (Color online) Schematic illustration of the pole structure 
and critical dynamics in the present model: a) the normal to super- 
radiance transition in the dissipative Dicke model, b) the normal to 
glass transition in the zero temperature equilibrium model, c) the nor- 
mal to glass transition in the dissipative model. 

a) When approaching the superradiance transition, two of the polari- 
tonic modes advance to the imaginary axis and become purely imag- 
inary before the transition point. This leads to the effective classical 
relaxational dynamics close to the transition. At the transition point, 
the Z 2 symmetry is broken by only a single mode approaching zero 
and becoming critical for J — > J c . 

b) For moderate disorder K, the poles are located on the real axis 
away from zero. For increasing K, the poles approach zero, with the 
closest pole scaling proportional to \K — K c \i, At K = K c the modes 
form a continuum which reaches zero and becomes quantum critical. 
No dissipative dynamics is involved. 

c) For moderate disorder K, the set of modes is located in the com- 
plex plane, away from zero. For increasing variance K, the modes 
get shifted closer to the origin, however, due to the scaling of real 
(oc \K - K c \ i ) and imaginary part (oc \K — K c \ 2 ), they neither become 
purely real nor purely imaginary. At K = K c a continuum of modes 
reaches zero. 



to < oj c , the atomic spectral density is described by 

Jl(co) - sgn(w) . 



2k(oJq+k 2 ^\oj\ 



(76) 



This unusual square root behavior of the spectral density in 
the glass phase, illustrated in Fig. [3] and also reflected in the 
pole structure Fig. 10 is a characteristic feature for glassy sys- 



tems that are coupled to an environment ifTTl [351 . It has been 
discussed previously in the context of metallic glasses, where 
collective charges couple to a bath of mobile electrons [ 1 1 1 or 
for spin glasses, where the spins couple to an external ohmic 
bath [35|. For intermediate frequencies, co > co c , the spectral 
density is linear, as it is known for the non-dissipative zero 
temperature case. In the limit k — > 0, co c is shifted to smaller 
and smaller frequencies, vanishing for k — 0. The cross-over 
frequency co c sets a time-scale t c — — , such that for times 
t < t c the system behaves as if it were isolated and one would 
observe the behavior of a T — quantum glass for (relative) 
time scales t < t c in experiments. On the other hand, the 
long time behavior, t > t c , of the atoms is described by over- 
damped dynamics, resulting from the coupling of the photons 
to a Markovian bath. This is a strong signature of low fre- 
quency equilibration of the atomic and photonic subsystem, 



which happens exactly at the glass transition (see Sec. V C 2 1 



C. Atom-photon thermalization 

We now discuss thermalization properties. The presence of 
quenched disorder in our model leads to an effective quartic 
atom-atom interaction, shown in Eq. ( |45j ), which allows for 
an energy redistribution to different frequency regimes. 



1. Atom distribution function 

In order to determine the atomic distribution function 
F., l (a>), we make use of the FDR (see Appendix IB] Eq. ( |B3[ >), 
which for the atoms described by a scalar degree of freedom 
simplifies to 



Q k (co) = FJco)(q r (oj)-Q a (oj)). 



(77) 



The atomic correlation function Q K is determined via 
Eq. ( |57] i. This equation contains the photonic Keldysh Green's 
function via A K (a>), and it is therefore evident, that the atomic 
distribution function will depend on the distribution function 
of the bare photons. The bare photon distribution function 
f ph (co) is again defined via the FDR, reading 

G*(w) = f ph {u) (G*M - G*(fii)) , (78) 

with the bare photon response and correlation functions 
GJ . Decomposing f-fs+fks into a symmetric 
/ s (w) = f s (-oj) and an anti-symmetric / AS (w) = -/ AS (-w) 
contribution allows us to rewrite A k (oj) in Eq. ( |57| i as 



2A*(w) = G$(a>) + Gf(-w) 

= 2/ AS (w) (AV) - AV)) + /.(to) (/(to) + g*(-oi) - g\aj) - g A {-aj)) 



CO 2 + k 2 + coi 
2|/ AS (to) + Vs(to) 



2ww, 



! o 



(aV)-a-V)). 



(79) 
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Inserting this result into the expression for the correlation 
function ( |57} , and making use of Eq. ( |58~} and its complex 
conjugate yields 



Q« e = (f^^f*)(Q R -Q A ) 



and thus identifies the atomic distribution function 



FJui) = / AS (w) + 






/ s (4 



(80) 



(81) 



This very general expression for the atomic distribution func- 
tion incorporates the two most important examples, either a 
coupling to a thermal or a Markovian bath. For the coupling 
to a heat bath, the bare photonic distribution function is fully 
anti-symmetric with f AS (oj) - cothf^j ,f s (u>) = 0, which im- 
plies that the atoms will be distributed according to a thermal 
distribution as well and experience the same temperature T 
as the photons. For the case of dissipative photons, the bare 
distribution function of the photons is fully symmetric, with 
f s (oj) - l,f AS (a>) = 0. Therefore the atomic distribution func- 
tion for this system is 



FJ.cS) 



CO- +K~ +COq 

2(jJll)Q 



(82) 



For small frequencies u> <K Jco^ + k 2 , the atomic distribu- 
tion function diverges as F. m (cl>) ~ -. This is the same asymp- 
totic low frequency behavior as for the thermal distribution 
function coth ( ^ J ~ — , such that for low frequencies, the 
system is effectively described by a thermal distribution with 
a non-zero low frequency effective temperature (LET) 



T efr = lim^. 



0)F a (co) 



4w 



(83) 



The atomic distribution F.„ and low-frequency effective 
temperature T dT in Eqs. ((82), ( |8"3") is identical to the distribu- 
tion function and LET of the photonic x-component, which is 
obtained by expressing the photonic action in terms of the x 
and p component, p - -j=(a' - a), and subsequently integrat- 
ing out the p component. This procedure is shown in App. [D] 
From the resulting action, the x component is described by a 
distribution function F xx (a>) = F M (u>), resulting from the cou- 
pling of the photons to the Markovian bath. Due to the strong 
atom-photon interaction, the atoms and the photonic x com- 
ponent equilibrate, resulting in the same distribution function 
and LET. 



2. Photon distribution function 

To compute the photon distribution function, we use the 
FDR 



G k {oj) = G R (u»F ph (u>) - F 9h G A {a», 



(84) 



which in this case is an equation for the 2x2 matrices G R ' A ' K 
and F, The matrix F solving Eq. ([84} is not diagonal, and the 
distribution of the excitations is determined by its eigenvalues 
f a . These are shown in Fig.|4]and illustrate the thermalization 
process of the system. In the normal and superradiant phase, 



the photons have a lower LET than the atoms, resulting from 
the frequency regime for which the dynamics is classical re- 
laxational. As for the spectral response, when the glass tran- 
sition is approached, this classical region is shifted towards 
a> - and the photon LET approaches the atomic effective 
temperature. At the transition, the photons and atoms have 
thermalized completely in the low frequency regime. 



D. Emergent photon glass 



In the glass phase, the condensate order parameter (a,) oc 
h Zz( cr f) = & vanishes for all photon modes (/)• However 
there exists a photon version of the Edwards-Anderson pa- 
rameter 

q EA = lim T ^ M ^ £fei<*;(* + t)x,(0> oc q EA , (85) 

where x = -4? (a + a') is the photon x operator and Eq. ( |85j ) 
only holds for the x-x correlations (and for those with fi- 
nite contributions to x-x). A non-vanishing photon Edwards- 
Anderson parameter implies an infinite correlation time for 
the photons, analogous to the atomic q EA -parameter. This is 
best illustrated by the correlator in the complex basis 

lim T ^ 00 <a(f + T)a ] (t)} = 1 lim T ^ 00 <x(f + -r)x(f)> = ^86) 

where we made use of the fact that the x-p and p-p correla- 
tions vanish for t — > co. Eq. ([86} implies that a photon which 
enters the cavity at time t has a non-vanishing probability to 
decay from the cavity at arbitrary time t + T, with r e [0, co). 
This highlights a connection to photon localization in disor- 
dered media 07] @8) . 

Close to the glass transition, the properties of the atomic 
system are completely mapped to the inverse photon Green's 
function. In the low frequency and small k limit, i.e. u, k «: 
(jJq,cd z , the inverse photon Green's function Eq. (64 1 has the 
expansion 



DIM = 



x-p 



2(oj 2 -A) 






-oj 



(87) 



such that the atomic low frequency physics is mapped to the 
photon x-x component. 

The determinant of D R _ vanishes at the zeros of D R , such 
that the photon propagator shows the same poles or branch 
cuts as the atomic propagator, and the scaling behavior at the 
glass transition obtained from the photons is identical to the 
one obtained from the atoms. The photon response properties 
induced by the atom-photon coupling are most pronounced in 
the x-x component G R X of the retarded photon Green's func- 
tion, 



G R Jco) = 



(io+ik) -\-2ti)QL R {(i))-ii)}. ' 



(88) 



For low frequencies, we can perform the same approximation 
as above to find 



G R x (co) = 



2{A-co 2 ) 






(89) 
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FIG. 11. (Color online) Spectral equilibration: Photon x-x spectral 
response 3K xx (<jS) = -2Im(Gf v (w)j in the glass phase for parame- 
ters K = 0.04, J = 0.12, id. = 2,co = 1,K = 0.02. As for the SR 
phase, it shows the same low frequency behavior as the atomic spec- 
tral response -2Im[Q R (a>)) (multiplied with a constant ^p). As for 
the atomic spectral density, one can clearly identify the overdamped 
regime with the unusual square-root behavior and the linear regime, 
separated by the frequency <±> c . 



Close to the glass transition and in the glass phase, the 
atomic retarded Green's function Q R determines the low fre- 
quency photon x-x response function. This is reflected in 
Fig. [TT] The discussion of the atomic response and scaling 
behavior in Sec. V B remains valid for the photons. 



E. Cavity glass microscope 

We now describe three experimental signatures (cavity out- 
put fluorescence spectrum, photon real-time correlation func- 
tion g (2) (r), and the photon response via driven homodyne de- 
tection) of the superradiant and glassy phases and their spec- 
tral properties. The cavity output is determined by the cavity 
input and the intra-cavity photons via the input-output relation 



„(a>) = y2K a(a>) + a.Jco), 



(90) 



with the cavity input annihilation operator ajco) and the aver- 
aged intra-cavity field 

m = r$ S l£ l a i {u>) (91) 

accounting for the M distinct cavity modes. 

1. Cavity output fluorescence spectrum 

The fluorescence spectrum S (u>) describes the (unnormal- 
ized) probability of measuring a photon of frequency to at the 
cavity output [33 1, and is defined by 



where a\ a (oS),a a JjjS) are creation, annihilation operators of 
the output field. Considering a vacuum input field, the fluo- 
rescence spectrum is expressed solely by the auto-correlation 
function of the intra-cavity field 

S(o) = <a t (w)a(w)) = / T e' w <a t (0)fl(T)> = /G < (w), (93) 

which is the "G-lesser" Green's function, occurring in the (±)- 
representation (see [61 68 1). Introducing also the "G-greater" 
Green's function 



iG > (oj) = Xc' c "<S(T)a t (0)>, 



(94) 



we can express response and correlation functions in terms of 
G </> 

G K (co) = G > (co) + G < (w) (95) 

G R (u>) - G A (co) = G > (co) - G < (co), (96) 

which yields 

s( w )^(gV)-cV) + g a m) 

= l - (GV) (F(oj) - 1) - (F(co) - 1) G A («)) ■ (97) 

In thermal equilibrium (F(a>) = 2«b(«) +1), where F is diag- 
onal in Nambu space, this expression simply reads 



S (oj) - n B {u))tR{oj) 



(98) 



and the fluorescence spectrum reveals information about the 
intra-cavity spectral density J[. 

In order to further analyze Eq. ( |97] >, we decompose the flu- 
orescence spectrum into a regular part and a singular part as it 
was done for the Keldysh Green's function in Eq. 



S(oj) = S„ g (o>) + 2nq EA 6(oj), 



(99) 



with the Edwards-Anderson parameter for the photons q EA . 
The regular part 5 teg is determined by the regular contribu- 
tions from the response and distribution function, G R ^ A ,F, 
which we have analyzed in the previous section. For small 
frequencies, F(oj) oc - and in the normal and SR phase, 
(G R - G A ) oc a), which leads to a finite contribution of S rcg to 
the spectrum. In contrast, in the QG phase, (G R - G A ) oc yfoj, 
such that 



StJoS) • 



l 



(100) 



has a square root divergence for small frequencies a> < co c QG 
(see Eq. |75|). This divergence is indicated in Fig.|7](c) and is 
a clear experimental signature of the glass phase. 

A further distinction between all three phases is possible by 
decomposing the fluorescence spectrum into a coherent and an 
incoherent part, where the coherent part describes the "classi- 
cal" solution (i.e. the part resulting from the presence of a 
photon condensate (5) + 0) and the incoherent part describes 
the fluctuations. Accordingly, the coherent part is 



(101) 



S(u)) = (ai,,(a>)a oul (w)), 



(92) 



S (OJ) = 27r|<a>| 2 (5(w) 
and the incoherent part reads 

S inc (cj) = S ms (cj) + 27r(q EA - |<a>| 2 )(5(w). (102) 
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Typical fluorescence spectra characterizing the three distinct 
phases are plotted in Fig. [7] For the normal phase, the spec- 
trum shows central and outer doublets associated with the hy- 
bridized atomic and photonic modes. Above the critical point 
for the superradiance transition, the doublets merge since a 
single mode becomes critical. However, compared to the 
single-mode transition, the central peak is much broader as 
a consequence of disorder. Additionally, in the superradi- 
ant phase, the fluorescence spectrum has a non-zero coherent 
contribution, which allows for a unique identification of this 
phase. 

In the glass phase, the doublets have merged after the emer- 
gence of a critical continuum of modes at 0) = 0, and one 
can clearly identify the square root divergence for small fre- 
quencies, as discussed above. Additionally, the singular be- 
havior of S(u>) in the glass phase is of incoherent nature, 
since (a"'") = 0. This combination of an incoherent zero fre- 
quency peak together with the absence of a coherent contri- 
bution uniquely defines the fluorescence spectrum in the glass 
phase and allows for a complete classification of the system's 
phases via fluorescence spectroscopy. 

The coherent contribution to the spectrum can be deter- 
mined via homodyne detection (see below), where (a) can be 
measured directly. 



2. Photon real-time correlation function g"\r) 

The time-resolved four-point correlation function of the 
output field 



g (2> (t,T) = 



( fl ou|M a ![( f +'>oul('+ T ) a oii|W> 

K«1('K U ,(0>I 



(103) 



reveals how the correlations in the cavity decay with the time 
difference t. In steady state, g (2) (t, t) only depends on the 
time difference t and we write g (2) (r). For t — > 0, g (2) (0) 
is a measure of the underlying photon statistics in the cavity, 
e.g. indicates bunching or anti -bunching of the cavity photons, 
respectively. 

In the open Dicke model, due to the effective temperature 
(cf. Fig. and Ref. ED), g (2) (0) > 1, describing photon 
bunching, as expected for thermal bosons. We find g (2 \Q) = 3 
for all the three phases, which stems from the off-diagonal 
atom-photon coupling in the Dicke model and coincides with 
the findings in Ref. [28 1 for the normal and superradiant 
phase. 

In the normal and superradiant phase, the long time behav- 
ior is governed by the classical low frequency dynamics, lead- 
ing to an exponential decay 

g (2) (r) ~ 1 + 2e- lKT . (104) 

This behavior is well known for the single mode Dicke model 
l28l and remains valid for the multimode case, away from 
the glass transition. In contrast, when the glass phase is ap- 
proached, the modes of the system form a branch cut in the 
complex plane and the correlation function in the glass phase 
decays algebraically, according to 



,(2) 



(T) ~ 1 + 



w 



where to = 0(1 /u>o). This algebraic decay of the correlation 
function provides clearcut evidence for a critical continuum 
of modes around zero frequency witnessing the glass phase. 
In Fig.p] we show g (2) (r) demonstrating this behavior. 

In order to compute the four-point correlation function 
( |103| l, we make use of Eq. ( |90"] i and the vacuum nature of the 
input field, i.e. the fact that all averages over a m , a in vanish. As 
a consequence, the operators for the output field in Eq. ( |103[ ) 
can be replaced by the operators for the averaged cavity field 
a, see Eq. ( [91) 1, The denominator in Eq. ( |103| > is then 



\{a\t)a(tW = \(a*Jt)a + (t)}\ 2 = IG^O)! 2 . 
The numerator similarly is expressed as 



(106) 



(a ] '(t)a ] '(t + r)a(t + r)a(t)) - (a*Jt)a*_(t + r)a + (t + t)2 + (i))7) 



Note that both expressions (Eqs. ( 106| l, ( |107| i) preserve the 
correct operator ordering of Eq. ( |103[ ), according to the differ- 
ent time-ordering on the (+), (-)-contour, respectively. 

The four-point function in Eq. ( |107| i can be expressed in 
terms of functional derivatives of the partition function Z. 
(Eq. ( |58| l) with respect to the source fields yt (Eq. ((37])). In 
the thermodynamic limit, the macroscopic action, Eq. ( |50"} , 
depends only on atomic and photonic two-point functions 
and, equivalent to Wick's theorem, the four-point function be- 
comes the sum over all possible products of two-point func- 
tions 

G {2 \t) = (a*_(t)a"_(t + T)a + (t + T)a + {t)) 

= {a*_(f)a + {t)) {a*_(t + T)a + (t + t)> 
+(a*_(t + T)a + (t)} (a*Jt)a + (t + t)> 
+{a_(t)a_(t + t)> <5 + (f + r)a + (f)> 
= IG^O)! 2 + |G < (x)| 2 + |G.:(r)| 2 , (108) 



with the anomalous G-lesser function 



G<(t) = -/<fl(r)a(0)>. 



(109) 



Inserting Eq. (jTOSJ into the expression for the four-point cor- 
relation function yields 



«®(t) = 1 + |^(T)| 2 + Pg 

with the two-point correlation function 

„(i),Vi - 2IM - (gtWggW 

•5 VJ G^O) (flt(0)8(0)>- 



(HO) 



(111) 



G < (t) is the Fourier transform of the fluorescence spectrum 
S (u>), as discussed in the previous section, and we therefore 
decompose it according to 



G^t) = q EA + G< k (t), 



(112) 



with G^ c X T ) being the Fourier transform of S^ico). In the infi- 
nite time limit, the regular part of G < (t) decays to zero, such 
that the infinite correlation time value becomes 



(105) 



a (X)( T \ — gEA _ 9EA 

r ^ M ?EA+C r *i g (0) t/EA+Nreg ' 



8 a \r) 



(113) 
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/V rcg denotes the occupation of the non-critical cavity modes. 
The way how this value of g (1) is reached in time is deter- 
mined by the -j= divergence of S lc „(w) for frequencies co < co c 
smaller than the cross-over frequency a> c , cf. Fig. 17] This 
leads to 



S (1) (T) 



^EA+Wrcg 



(?)'■ 



(114) 



where r c = — and f o has to be determined numerically. This 
algebraic decay to the infinite t value of the correlation func- 
tion with the exponent v - \ has also been found in Ref. ||34| 
for the correlation function of a spin glass coupling to a finite 
temperature ohmic bath, in line with the discussion of uni- 



versality in Sec. V B The finite temperature exponent results 



from the non-zero effective temperature of the system, which 
influences the correlation function. For the case of T dY - 
this exponent changes to v = | but the spectral properties are 
left unchanged. 

The non-zero value of the two-point correlation g (1) (Y) — » 
_ ^ A N f or t — > oo serves as a possible measure of the pho- 
tonic Edwards-Anderson parameter q EA in the glass phase: q EA 
can be inferred from a correlation measurement, if the total 
photon number in the cavity N m = q EA + N m has been mea- 
sured separately. 

Taking the absolute value of g (1) (r) in Eq. 
the dominant contribution 



FBI), leads to 



I^W 



\qEA+N, ce ) 



+ 2 



ggA 



C1EA+N, 



-M 



(115) 



as displayed in the asymptotic behavior of the four-point cor- 
relation function (Eq. ( |105| l , where we have absorbed the 
prefactors in the definition of tq and normalized the long-time 
limit to unity. 

While the non-zero value of g (1) (r — > oo) is caused by crit- 
ical poles of the system, it does not include any more infor- 
mation about the pole structure of the system and may for 
instance be caused by a single critical pole, as it is the case for 
the superradiance transition. However, the algebraic decay to 
the infinite correlation time value of g (1) , and the same for g (2) , 
is a clear signature of a branch cut in the complex plane and 
therefore a continuum of modes reaching to zero frequency. 
This in turn is a strong signature of the critical glass phase in 
the cavity. 



3. Photon response via driven homodyne detection 

Here we relate homodyne detection measurements of the 
output signal to the quadrature response functions in the 
Keldysh formalism and calculate the corresponding signal. 
This gives predictions for the experimental analysis of the 
spectral properties and the scaling at the glass transition, 
which have been discussed in previous sections. 

In the process of homodyne detection, the output field a M is 
sent to a beam-splitter, where it is superimposed with a coher- 
ent light field P(t) - p e -«H8<+fl) w j t } 1 frequency cop, amplitude 
/3 and phase 9. After passing the beam-splitter, the intensity 
of the two resulting light fields is measured and the difference 
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FIG. 12. (Color online) Illustration of homodyne detection of a 
weakly driven cavity. The cavity is driven via a weak coherent input 
field 77(f) entering the cavity through one of the mirrors. Then a ho- 
modyne measurement is performed on the output signal of the driven 
cavity. For this, the output signal is superimposed with a reference 
laser /3(f) via a 50/50 beam-splitter and the difference current of the 
two outgoing channels is measured. From this, the response function 
of the photons in the cavity can be measured by tuning the relative 
phases and frequencies of /3(f) and 77(f), as explained in the text. 



in this measurement (the difference current) for the case of a 
50/50 beam-splitter is described by 

n-(t) = i(ai t (t)/3(t)-0*(t)a Mt (t)) 



).(H6) 



Here, we added a conventional phase shift <f> — | of the beam- 
splitter. For the case of a vacuum input field, Eq. ( |116| > sim- 
ply measures the steady state expectation value of the cavity 
quadrature components 



X $ * „,(*) = e^mty^ + e 



->'(<K);}f 



cv{t)e 



-WJjit 



(117) 



with the intra-cavity operators 5, as defined in Eqs. ( |90"| i, ( |9Tj >. 
This quantity indicates a finite superradiance order parameter 
(5), but for the steady state contains no further information. 

This situation changes when the input field is changed from 
the vacuum state to a weak coherent laser field 77(f). For this 
special case, the difference current in Eq. \\\6\ is modified 
according to 

n.(i) = i (77*73 - fTrfi (f) + i ^FTk (a\t)P(t) - p*(t)a(t)) 

= i (77*73 - 73*77) (f) + V243| (X e _|^(0) ■ (118) 

For the special case of the input field coming from the same 
signal as the reference laser, we have 77(f) = fi(t) (which 
we assume from now on for simplicity) and the first term 
in Eq. ( |118| l vanishes. The main difference here, is that the 
quadrature operator Xg M/1 (t) is not evaluated for the steady 
state but for a state which has been perturbed by the weak laser 
field 73(f). For a weak laser amplitude |/3| <K 1, the system stays 
in the linear response regime and the difference current is pro- 
portional to the retarded Green's function for the quadrature 
component Xg + * tW/J as we proceed to show. 
The interaction between the cavity photons and the radiation 
field outside the cavity is commonly described by the Hamil- 
tonian 



H m — i y2i< I a' a m — a m a\ , 



(119) 
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which after a transformation to the Keldysh action and replac- 
ing the input fields by the coherent light field j3(t) enters the 
action as 

S m = VSJk(a;,a;)(a))Hr*(| c )(to) + hx., (120) 

which is exactly the form of a source term in quantum field 
theory, generating all Green's functions of the system via 
functional derivatives with respect to the fields /3. Express- 
ing the action ( |120| > in terms of Keldysh components of the 
quadrature fields Xg^ yields 

S in , = ^TKl(x cfi -^,X qfi -^)(cS)(r x l ^| Wl21) 

The linear response (X fl _» iW/) ) (1) (f) of the quadrature expecta- 
tion value is expressed as (see appendix |C|) 

(X e -^) m (t) = -2kW X G£_ (* - O, (122) 
with the quadrature response function 

G x 8 _ ?; ^( f " O = -»W " f){[x e - h a > ,(t),X e . l ^(0]lX23) 
For the specific choice of 6 — |, 

X^ 5ju . = V2x„„ = (g(0«**' + ~a\t)e-^*) , (124) 



.Wy3 



"1»|] 



the response function Gv (t-t')-2Gl (t - f) becomes 

the x-x retarded Green's function in a frame rotating with the 
laser frequency cop. In this case, the difference current is a 
direct measurement of the x-x response 



n-« = -443| 2 £Gl(r-0, 



(125) 



which we have discussed in detail in Sec. V D The frequency 
dependence of x Wfl , indicated by the subscript o)p, coming 



from Eq. ( jll7[ ), can be used to scan through different fre- 
quency regimes and directly access the atom and photon x-x 
spectral response. 



VI. CONCLUSION 

We have developed the non-equilibrium theory of the mul- 
timode Dicke model with quenched disorder and Markovian 
dissipation, and provided a comprehensive characterization of 
the resulting phases in terms of standard experimental observ- 
ables. The main theoretical findings relate to the interplay 
of disorder and dissipation. We establish the robustness of a 
disorder induced glass in the presence of Markovian dissipa- 
tion. This concerns, for example, the presence of an Edwards- 
Anderson order parameter and the algebraic decay of corre- 
lation functions in the entire glass phase. Central quantita- 
tive aspects, such as the decay exponents of the correlation 
functions, are strongly affected by the presence of dissipation. 
Disorder leads to enhanced equilibration of the atomic and 
photonic subsystems for both the spectral (response) and their 
statistical properties. The spin glass physics of the atoms is 



mirrored onto the photonic degrees of freedom. We presented 
direct experimental signatures for the atomic and photonic dy- 
namics that allow unambiguous characterization of the vari- 
ous superradiant and glassy phases. 

Several directions for future work emerge from these re- 
sults. In particular, the realization of disorder may not be gov- 
erned by an ideal single Gaussian probability distribution in 
experimental realizations of multimode Dicke models. This 
may concern, for example, effects relating to the finite num- 
ber of cavity modes (M) or effective two-level atoms (N). 1 /N 
corrections contain information on the critical behavior close 
to the conventional Dicke transition J27l[30l . with similar fea- 
tures expected for the glass transition. While we expect the 
main glassy features to be robust to such finite-size effects, it 
would be interesting to study a concrete cavity geometry with 
specific information of the cavity mode functions. 

Furthermore, with our focus on the stationary state we did 
not touch upon the interesting questions of glassy dynamics 
[35, 63 1 in this work (for thermalization dynamics of the sin- 
gle mode Dicke model, see (69)). An interesting problem is a 
quantum quench of the open, disordered system. In particular, 
the non-universal short time and transient regimes should con- 
tain more system specific and non-equilibrium information. In 
the long time limit, the nature of aging and dependencies on 
the aging protocol remains to be explored. 
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Appendix A: Photon fields for superradiant phase 

In order to describe a system where the particle number is 
not conserved, as it is the case for the photons in the Dicke 
model, we introduce the spinor field 



A (t) - I a "J {t) 



(Al) 



containing the bosonic fields a a .(f), a* a ft) for a quantum state 
j and with index a - q,c. The corresponding adjoint field is 



A 1 " (f) = (a* (f),a .(?))■ 



(A2) 
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The action for a quadratic problem is (for simplicy we con- 
sider only a single quantum state) 



S = f lt ,(AUt),Al(t))D 4x4 (t,0 



AJf) 



(A3) 



where 







^^'''-l^oKnl- 1 ^ 11 ''-''^ 



is the inverse Green's function. The Keldysh correlation and 
retarded Green's function are also 2x2 matrices, which can 
be expressed in terms of operator averages according to 



G R 2x2 (t,o = (D R 2x2 y\t,t') 



(A5) 



and 



■ft(t-t'M\ WWCO] [a(t),a(t')] 
Ml [a\t),aHt')] [aHt),a{t')] 



G K lxl (t,f) = - (G R 2x2 o Df x2 o Gj x2 ) (t, t') 

{a(t),aUt')} {a(t),a(t')} 
{aHt),a ] '(t')} {a\t),a(t')} 



(A6) 



In Eq. ( |A6[ ), the o-operation represents convolution with re- 
spect to time. 

For the Dicke model with strong atom-photon coupling, it 
is reasonable to transform to the x-p representation in terms 
of real fields 

x(t) = Jg (a*(t) + a(t)) , p(t) = ^ (a*(t) - a(t)) . (A7) 
This is done via the unitary transformation for the fields 



(x(f),p(t)) = (a?(t),a(t))--l j ! 



V2 



(A8) 



and the Green's function 



The same can be done for the advanced and Keldysh Green's 
functions, leading to the expressions for response and correla- 
tion functions as discussed in the main text. 



Appendix B: Markovian dissipation vs. quenched disorder 

As anticipated in the main text, the quenched bath, resulting from the coupling to a static distribution, is fundamentally 
different from the Markovian bath, represented by the fast electromagnetic field outside the cavity. While the dynamics of the 
quenched bath is frozen on the time scales of the system, the dynamics of the Markovian bath happens on much faster time scales 
than those of the system. As we will see, both types of bath inherently lead to non-equilibrium dynamics of the system since the 
system-bath equilibration time becomes infinite. For both cases this implies a non-equilibrium fluctuation-dissipation-relation 
(FDR), connecting response and correlations via a non-thermal distribution function. 



1. Non-equilibrium fluctuation dissipation relation 



Correlation and response properties are not fully independent of each other but connected via fluctuation-dissipation relations, 
which we will briefly introduce in this part. In a system with multiple degrees of freedom, the response properties are encoded 
in the retarded (advanced) Green's function G R( - A \t, f), which is defined as 

G R (t,f) = -ie(t-t){[a i (f),a)(f)]), 



(Bl) 



with the commutator [•,•], the system creation and annihilation operators a\,a j and G R (t,t') - \G A (t',ij\ . The correlation 
function on the other hand 



C ij (t,f) = {{a i (t),a)(f)}) = iG*(t,f) 

is defined via the anti -commutator {■, ■} and defines the Keldysh Green's function G K (t, t') 
The fluctuation dissipation relation states 

G k (oj) = G R (aj)F(a>) - F{aJ)G A (a>) 



(B2) 



(B3) 



and relates the response and correlations of the system via the distribution function F(a>). In thermal equilibrium, the distribution 
function is fully determined by the quantum statistics of the particles and the temperature T according to 



Fij((o) = 6ij(2n B (co)+ 1), 



(B4) 



with the Bose distribution function n B . As a result, in equilibrium, it is sufficient to determine either response or correlation 
properties in order to gain information on each of these. 
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2. Effective system-only action 

In this part, we present a derivation of a system-only action after elimination of the bath variables via Gaussian integration. 
Depending on the nature of the bath, different distribution functions will be imprinted to the system. We start with the general 
action of the bath, which we consider to be well described by a quadratic action and in the (±)-basis 

s B = E, {,, {dM M) ( % %- ) ' ('. O ( f_^] ) . (B5) 

with the bath variables ^ and the bath mode index //, which will be chosen a continuous index below. The Green's functions for 
the uncoupled bath variables are assumed to be in equilibrium and read 

G;-(t, f) = G<(t, t') = -jn(^) e -''^ ( '- f,) (B6) 

G~\t, t') = G>(f, t') = -i{n{u,) + 1) e-Hfr-O (B7 ) 

G; + (t, t') = Gfa t') = 0(t - f)G> + 0(t' - i)G< (B8) 

G"(t, t') = G*(f, t') = 9(t - f)G< + 9{t' - i)G>, (B9) 

with the bath frequencies co^ and the familiar Green's functions G-lesser, G-greater, the time-ordered and the anti-time ordered 
Green's function. The linear coupling between system and bath is 

5, = Z„ <fn 5,(al<f)Mi))[ I _°! if f^ ) + h.c, (BIO) 

where a' ,a are the system's creation and annihilation operators. For simplicity we consider only a single quantum state of the 
system, but a generalization to many states is straightforward. The partition function is of the general form 

= J D[a,a ] ]e^ U D[^,^]e iiS ' +s A , (Bll) 

e' 5 eff 

where S s is the bare action of the system. Now we integrate out the bath via completion of the square. The contribution iS e sji 
of the //th mode to the effective action reads 

c r fi CfUA Utwl GfM G^COW « + (0 \ mi9 , 

S^[a,a^=7 tl JJ<(f),-al(t))^ ^_ ((>(/) G ^^ n j[_ a(f) J- (B12) 

The signs for the operators on the - contour come from the backward integration in time. Thus the mixed terms will occur with 
an overall - sign, while the ++ and — terms come with an overall +. Summing over all the modes /i we obtain the effective 
action of the bath for the field variables of the subsystem. We now take the continuum limit of densely lying bath modes, centered 
around some central frequency «o and with bandwidth §. That is, we substitute the sum over the modes with an integral in the 
energy Q weighted by a (phenomenologically introduced) density of states (DOS) v(Q) of the bath 

Z» y, - IT/ dn r (0)v(ii) (b i 3) 

and obtain 

S^a, at] = - JT* XWWm £(4(0, -aim ( gig J|$ ) ( ^ } ) , (B14) 

where in addition we have used the translation invariance of the bath Green's function, G^it, f) - G°^(t-t') to suitably shift the 
integration variables. Eq. ( |B14[ ) is a general expression for an effective system action resulting from a coupling of the system to 
a bath of harmonic oscillators with a coupling that is linear in the bath operators. In the case of a strong separation of time scales, 
the effective action can be further simplified. Here we consider two extreme and opposite limiting cases, namely a Markov and 
a quenched disorder bath. 
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3. The Markov approximation 

The Markov approximation is appropriate when there exists a rotating frame in which the evolution of the system is slow 
compared to the scales of the bath, i.e. o> sy! <K uq, §, such that the system is considered as being static on the typical time scale 
of the bath. This leads to a temporally local form of the resulting effective action. As an example, we derive the (±)-part of the 
effective action 



Eq. |B6| . 

- iyv 



dtaUol drj — y(Q)v(Ci)G^(j)a + (t - r) 

-yv I dtal(f)i J dr j ■^G^(T)\a + (t- S ) 

f dt al{t) Ifdrf — «(Q)e-' (n -"» )T ] a + (t. s ) 



ss Tim I dt aL(t)a + (t-s). (B15) 

In the second line, we made use of the Markov approximation, i.e. the time evolution of the system is much slower than the one 
of the bath in the rotating frame, and the coupling and DOS are constant over the relevant frequency interval. In the third line, we 
replaced the Green's function by its definition (in the rotating frame). Finally, in the last line, we introduced the particle number 
h - n(w ( )) at the rotating frequency and the effective coupling 2k - yv. Performing these steps for all the four contributions to 
the action in the (±)-basis leads to the action 

SJLa,eP] = fdt(at(t),ai(t))^( a a + ® J, (B16) 

which is local in time, containing the Markovian dissipative self-energies 

. I 2h+ 1 -2(n+ 1)\ _,,_ 

^ = ' K [-2n 2n + lj- < B17 > 

Transforming this self-energy to the Keldysh representation, we finally obtain 

^ = i«(H l4n \ 2 )- (B18) 

The additional contribution to the distribution function F Mar (w) for the Markovian case is obtained from the FDR for the self- 
energies 

T, K (tu) = F(oj) (E R (w) - S-V)) . (B19) 

For the case when the system couples only to the Markovian bath or to an additional thermal bath, these contributions are 
infinitesimal and only those from the Markovian bath have to be taken into account, yielding 

iK(4h + 2) = F{u)2i K , (B20) 

i.e. the distribution function 

F(oS) = 2n + l. (B21) 

In this expression, the frequency dependent particle distribution n(u>) has been replaced by the relevant particle number n(a>o) 
of the bath. The interpretation of this, is that the dynamics in the bath are so fast compared to the system, that the for the full 
frequency regime, the system only couples to the slowest bath modes (in the rotating frame), located at u> - coc*. This makes it 
impossible for the system to equilibrate with the bath and it can therefore not be described by a thermal distribution, i.e. stays 
out of equilibrium. 

4. The quenched bath 

The quenched bath is located in the opposite limit of the Markovian bath, i.e. it constitutes of a system bath coupling, such 
that there exists a rotating frame for which the system dynamics is much faster than the bath dynamics, i.e. wq, ■« w, ys . The 
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corresponding approximation is to assume that the bath is static on the relevant time scale of the system and the resulting effective 
action for the system is infinite range in time. In this case, the contribution to the action for the (h — )-component reads 

dtaliOj drj — y(D.)v(Ci)G^(T)a + (t-T) 

dtj dral(t)\\ -^n( 



lyv 



(Q) \a + (t-T) 



' da> 



HkN I — a' (a>) d(co) aJaS). 

2n 



(B22) 



In the second line, we inserted the definition of the Green's function and made the approximation of a slowly varying bath as 
well as a constant DOS and coupling, v, y. In the third line, we replaced yv - 2k and inserted the average particle number of the 
bathJV. 

Repeating these steps for all contributions to the action in the (±) basis and subsequently transforming to the Keldysh basis, 
we have the self-energy 



E Q (w) = iKd((D) 







1 



-1 4N + 2 • 



(B23) 



This contribution is structurally different from the one from integrating out the Markovian bath, since it only acts at a> - 0. As a 
result, the distribution function for the system is only changed for a> - compared to the uncoupled, bare system. And therefore 



F(a>) = 



2N+1 if 0) = 
F baie (w) ifw*0 ' 



(B24) 



where F bm is the distribution of the bare system. In contrast to the Markovian case, where we obtain a constant distribution 
for all frequencies and therefore higher system frequencies are strongly pronounced, the quenched bath shifts the occupation 
distribution to the very slowest modes of the system, therefore implying very slow dynamics on the system. This is reflected in 
the modified FDR and the appearance of a glassy phase, as discussed in Sec. |IVB] 

The picture obtained from these extreme cases of possible system bath couplings is quite transparent. For an equilibrium 
system, one assumes that the bath is such that for any possible frequency of the system, there exists a continuum of modes in 
the bath, such that thermalization of the system will happen on the whole frequency interval. In contrast, when the bath modes 
are located at much higher frequencies than the system, all the system modes interact the strongest with the slowest bath modes, 
leading to a distribution function as depicted in Eq. ( |B21[ ) and avoiding direct thermalization. On the other hand, for a bath that 
evolves on much slower time scales than the system, the picture is reversed, and only the slowest modes of the system interact 
with all the bath modes in an equivalent way. For the extreme case of a static bath, all the bath modes interact with the system's 
zero frequency mode, and the distribution function becomes the one in Eq. ( |B24[ ). This is again a non-equilibrium distribution, 
such that the system does not directly thermalize. 



Appendix C: Linear response in the Keldysh formalism 

A common experimental procedure to probe a physical sys- 
tem is to apply a small external perturbation and measure the 
system's corresponding response. If the perturbation is suffi- 
ciently weak, the measured response will be linear in the gen- 
eralized perturbing force. Here we review this construction 
in the Keldysh formalism in order to provide the background 
for the connection to the input-output formalism of quantum 
optics made in the text. 

We consider a setup, where the hermitian operator O = & 
is measured after a perturbation of the form 



HJt) = F(f)0 



(CI) 



has been switched on at t — 0. Here, the (unknown) real val- 
ued field F(t) oc 0(f) is the corresponding generalized force. 
The expectation value 



(6)(t) = |Tr [p{t) 6) 



(C2) 



is evaluated by introducing a source field h(t), such that 

« 5 >W = iSL. (C3) 

where 

Z(/j)=Tr(e-^ + /*^)<5®). (C4) 

Expressing Z in a real-time Keldysh framework, we have 

Z(h) = JD[if/*, </r] e ''W-« e iSS\h,rm t (C5) 

where So is the unperturbed action and {if/, if/*} are the complex 
fields representing the creation and annihilation operators of 
the system (in the ±-basis). The term 

6S[h,if/*,if/] = fdt(h+(t)0+(f)[f + ,i/f+] 

-h-(f)0-(f)W-,f-\) (C6) 

contains the source fields h ± coupling to ± which polynomi- 
als in if/*, if/. The expectation value ( |C2[ ) transforms according 



to 



(0(f)) = <O + (0> = <O_(0> 



i<O + (0 + O.(0>, (C7) 



whereas the averages on the right always mean averages with 
respect to the functional integral. In terms of functional 
derivatives of the partition function, we find 



(0(f)) = - 



2 \5h+(f) 6h-(t) 

i 5 



Z(h) 



/?=o 



V2 Sh q (f) 



Z(h) 



(C8) 



/)=o 



The second equality results from a rotation to the RAK rep- 
resentation and determines the time-dependent expectation 
value of 0(t) for a system described by the action Sq, In order 
to incorporate the perturbation ( |Cl| l, we add the perturbation 
to the bare action of Eq. ( |C5| l 
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Now we can expand the expectation value of to various 
orders in the force. The zeroth order simply is the expectation 
value in the absence of the perturbation: 



0(t))<® = - -L jX- Z(h ,F)\ 

S2M q (t) v '\F=h=0 



The linear order term is then obtained via 



(CIO) 



(0(f)) 



6F + (f) ) F=Q 



1} = J dt'F + (f)i 

+F -(O(^^«5(0>) . (Ci 

\SF-(f) ) Fsd0 



1) 



So — > So + J dt (F + (f)0 + (f) — F.(f)0-(f)) . (C9) which after a translation into the RAK representation reads 



(0)| +F.(f) 

F=0 



--U dt ' F( 4^ + ^) {o+(t)+o - (t) L 



— — (0 + (t) + O_(0> 
6F.(f) /F=0 



I J \SF + (f 

r , , s 2 

dt'F(f) 

J SF c (t')5hJf) 



) SF.(f) \Sh+(t) Sh.(f) 



Z(h, F) 



Z(h, F) 



F=h=0 



h=F=0 

jdt'F(t')G R 00 (t,t'), 



(C12) 



where we made use of (at the point where we extract physical 
information) F+(f) - F-(f) = F(f), and furthermore that f < t, 
such that the last equality indeed yields the retarded Green's 
function for the operator O. The integral in ( |C12| i runs from 
minus infinity to plus infinity, whereas the retarded Green's 
function defines the upper bound being t and the force F(f) 
sets the lower bound to be to since it vanishes for t < to when 
the perturbation is switched on at t — to- Since the integral for- 
mally runs from minus infinity to plus infinity, we can switch 
to frequency space, where for the time-translational system 
(stationary state) we find 



«5>< V) = -F(o>)G R nn (aS) 



(C13) 



1. Example: laser field induced polarization of cavity atoms 

The polarization of an atomic two-level system can be ex- 
pressed as 



P(t) = (fi R cr*(t) + mcry(t)), 
or after a rotation around the z-axis 

P(f)= l x((T x (f)). 



(C14) 



(C15) 



r 



We are interested in the response of the polarization to a per- 
turbation of the system by a coherent monochromatic light 
field. Since the coupling of the light field is proportional to 
the polarization, the corresponding Hamiltonian reads 



H(f) = Q.(f)a- X , 



(C16) 



where Q(f) = 6(t)yiE(t) is the generalized force and E(t) is 
the electric field. The corresponding action for this problem 
is then 

S =S + SS[h,£l,(f>] 



So+ I dt h q (f)<f> c (t) + h c (t)<f> q (t) 
+Q q (f)<p c (f) + &c(t)<P q (f), 



(C17) 



where we have replaced <x v by the real fields cf> as in Sec. IV A 
Applying ( |C12[ ), we then find 

P , , s 2 

<H I dt£l(t 



P w (t) 



6Q c (f)Sh q (t) 

= -n | dt'D.(t')Q R (t - f), 



Z(h, O) 



n=/i=o 



(C18) 



where Q R (t - f) is the retarded atomic propagator as in the 
previous sections. Now we again switch to frequency space 
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and use the definition of Q, such that we find 
P (l) (io) = /?£(w)G*(w), 



(C19) 



where we have absorbed the 0-function into the electric field. 
This equation identifies the retarded atomic Green's function 
that we used in the previous section with the linear atomic 
susceptibility ?c\(o), which is commonly used in a quantum 
optics context. 



Appendix D: Distribution function of the photon x-component 

In this section, we derive the distribution function for the 
photonic x-component and show that it is identical to the 
atomic distribution function, proving that the atoms equili- 
brate with the photon x-component. 

The Keldysh action describing the bare photon degrees of 
freedom is given by Eq. pi} and we express this action di- 
rectly in the Nambu basis, using the vector 



A,(w) 



' a c (oj) ^ 
a* c (-u>) 
a q {co) 
a*(-w) ) 



(Dl) 



the photonic action reads 

S ph = £ A+ (w)D 4x 4(w)A 4 (w), (D2) 

with the inverse Green's function in Nambu representation 



D 4x4 (w) 



2 



x2 



(co + iK)cr z + wo 1 2x2, 



{(jJ-iK)<T z + W 1 2x2 



2M 



2x2 



(IP3) 



The action ( |D2| i can also be expressed in terms of real fields 
by performing the unitary transformation 



x a {co) \ _ j_ I 1 1 a a {cj) 
Pa(oj)j ^2\ i ~i)[a* a (-to) 



(D4) 



with a - c,q. After this transformation, we express the action 
in terms of the real field 



V 4 (cj) = 



( X c (LO) \ 
Pc{m) 

x ? (w) 

Pq(u) ) 



such that 

S ^ = L VZ(-(o)D x . p (co)V 4 (u), 
with the inverse Green's function 



D x - p (oj) = 









-w 


k — ico 








—k + ico 


-w 


-Wo 


—K — icO 


Hk 





+ ico 


-w 





2'lK 



(D5) 



(D6) 



(D7) 



The action ( |D6| l is quadratic in the fields x a and p a and we 
can eliminate the p-fields from the action via Gaussian inte- 
gration. The resulting action is 



with the field 



s-il f (*(*)- 



X(a>) = 



x e (w) 
x„(w) 



and the inverse Green's function 




D x {co) 



(w - we) 



(W + iK) — Wy 

bj 



2 \ 



(D8) 
(D9) 

(D10) 



The distribution function F x (co) for the x-field is obtained via 
the fluctuation-dissipation relation 



D K X (w) = F x {co) (D*(w) - DiicS)) , (Dll) 



yielding 



FJco) 



llOQiO 



(D12) 



This is indeed identical to the atomic distribution function, 



that we have computed in Sec. VC which proves that the 
atoms equilibrate with the photon x-field. 
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